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THE DEVELOPMENT OF ACCURATE AND EFFICIENT 
METHODS OF NUMERICAL QUADRATURE* 

By 

T. Feagin 
ABSTRACT 

Some new methods for performing the numerical quadrature of an 
integrable function over a finite interval are described. Each method provides 
a sequence of approximations of increasing order to the value of the integral. 
Each approximation makes use of all previously computed values of the inte- 
grand. The points at which new values of the integrand are computed are se- 
lected in such a way that the order of the approximation is maximized. The 
methods are compared with the quadrature methods of Clenshaw and Curtis, 
Gauss, Patterson, and Romberg using several examples. 
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THE DEVELOPMENT OF ACCURATE AND EFFICIENT 
METHODS OF NUMERICAL QUADRATURE* 


INTRODUCTION 

One of the most basic problems in numerical mathematics concerns the 
numerical computation of the definite integral of an integrable function defined 
over some finite interval. Most numerical methods that have been devised for 
performing such quadratures approximate the value of the integral by a linear 
combination of function values, that is, 


I = 


f f (x) dx = £ W i f(x.) 
i=l 


where f(x) represents the function defined on the interval (a, b); w A denotes the 
weight of the value of the function at the i th point or abscissa, x.; and E n 
designates the error committed when the approximation is made. 

There have generally been two basic approaches to the problem of numeri- 
cal quadratures, adaptive and non-adaptive. If the abscissae are selected at 
some intermediate stage of the process depending upon the behavior of the 
integrand, then the pi cess is said to be adaptive. If the process evolves inde- 
pendently of the behavior of the integrand, the process is said to be non-adaptive. 

*This research was conducted while the author held a NRC Resident Research Associateship at 
Goddard Space Flight Csr.i jr, Greenbelt, Maryland. 
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The methods developed here may be used readily in either approach, because an 
estimate of the error is available at each stage. The error estimate could be 
used to make a decision regarding the subdivision of the interval into smaller 
parts and the formulae could then be used on the smaller intervals. 

If the behavior of the integrand is sufficiently well understood at the outset, 
it may be feasible to determine a value of n for which the error E n is accep- 
table. In such a context, the Gaussian formulae generally provide the most ac- 
curate and economical results. The Gaussian weights and abscissae are 
tabulated, for example, by Stroud and Secrest (1966). 

On the other hand, if little is known about the behavior of the integrand, it 
may be necessary to compute many approximate values of the integral before 
an acceptable result is found. Furthermore, it is usually necessary to determine 
whether or not a given result is acceptable. It is therefore recommended that 
any practical procedure for performing quadratures be able to provide a sequence 
of results of increasing accuracy, providing simultaneously reliable error esti- 
mates which may be used to decide when to terminate the procedure. The 
methods presented here have the above characteristics and also appear to be 
very efficient with regard to the number of integrand values required in order 
to achieve a given accuracy. 

One method of numerical quadratures that has the ability to provide a 
sequence of approximations of increasing order would be to use a sequence of 
Gaussian formulae and to estimate the error as the difference between the 
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is evaluated include all old points, thereby causing the algorithm to be quite 
economical. It should be remarked that the additional effort required in order 
to obtain the coefficients of the series is only necessary if the error estimate 
or an indefinite integral is desired. 

It is important to realize that no estimate of the error E n is completely 
foolproof. There is always some function which will cause any estimate to be a 
gross underestimate or a gross overestimate of the true error. Usually the 
effort expended in order to improve the estimate of E is equivalent to improv- 

n 

ing the approximate value of the integral itself. The most essential feature of 
a good quadrature method is therefore the ability to obtain a sequence of accu- 
rate results with as few function evaluations as is possible. 

Patterson (1968) has derived a set of quadrature formulae which have 
several desirable features. The method is comprised of a sequence of quad- 
rature formulae (beginning with the midpoint rule, the Gauss three-point rule, 
and Kronrod's seven-point rule) each of which contains all previous abscissae 
and the new abscissae and the weights of each formulae are selected so as to 
optimize the degree of the formula. The new set of n-point formulae are of 
degree (3n + l)/2. These formulae may be differenced to provide an estimate of 
the error, or perhaps a series of Chebyshev polynomials could be fitted to the 
values of the integrand and integrated term-by-term and the last few terms of 
the series could be used as an estimate of the error. 
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results obtained for the different orders. The drawback with such a method is 
that the Gaussian formulae have few abscissae in common for different values 
of n, thereby causing the total number of function evaluations required to obtain 
an acceptable result to grow rapidly with the order of the approximation. 

Another method which might be used to generate a sequence of results of 
increasing order is based upon the use of Newton-Cotes formulae. In this case, 
the formulae do have many abscissae in common for different values of n. The 
problem here is that the higher-order Newton-Cotes have undesirable features 
with regard to round-off errors and convergence properties (Hildebrand, 1956). 
Romberg (1955) has developed a method which is also based upon equally- spaced 
points that avoids these difficulties. However, the Romberg method achieves 
high-order results very slowly as the number of function evaluations increases. 
Bauer, Rutishauser, and Stiefel (1963) have succeeded in reducing somewhat the 
number of function values required by modifying the Romberg algorithm, but the 
modification causes the method to be more susceptible to round-off errors. 

Clenshaw and Curtis (1960) recommend a method of numerical quadrature 
based upon the development of the integrand in a Chebyshev series of polynomials. 
The quadrature is accomplished by integrating the series term-by-term. The 
error is estimated by examining the last few terms of the integrated series. If 
the coefficients of the series are converging rapidly enough, the last few terms 
generally provide a conservative estimate of the error. This method is 
particularly advantageous if an indefinite integral is desired. If the number of 
terms in the series is increased from N to 2N-1, the points at which the function 
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DESCRIPTION AND DERIVATION OF THE METHODS 


The methods presented here are founded upon the same basic philosophy as 
that adopted by Patterson (1968). Such methods are characterized by the initial 
formula or rule to which all subsequent abscissae are added and by the way in 
which these points are added. In the same way that Patterson's method has the 
midpoint rule as its foundation, the two methods developed here are based upon 
the trapezoid rule and the Gauss two-point rule respectively. 

The first method is based upon the trapezoid rule and is comprised of a 
sequence of quadrature formulae with 2+1 abscissae for k = 0, 1, 2, . . . 

The weights and abscissae of the first five members of the sequence are given 
in Table I, where the interval (a, b) has been chosen to be (-1, 1) for simplifica- 
tion. The error term of each formulae of the sequence is also tabulated. 

The second method is based upon the Gauss two-point rule and is com- 
prised of a sequence of quadrature formulae with 3 • 2 k -1 abscissae for 
k = 0, 1, 2, . . . . The weights, abscissae, and the error term for each of the 
first four members of this sequence are given in Table II. 

The formulae may be derived in the following way. For each formulae of 
a given method, the n abscissae may be characterized as the roots of an n th 
degree polynomial, G n (x), where 

n 

G n( X > = V j - 

j=0 
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Table I 


Method I 


Two-point Formula of Degree One 

(Trapezoid Rule) 

E n = -6.67 x 10 _1 f (2 te) 

± X. 
1 

w 

1 



1.0000 0000 0000 

1.0000 0000 0000 



Three-point Formula of Degree Three 

(Simpson's Rule) 

E =-1.11 

n 

x 10- 2 f (4 \f) 

± x. 
1 

w. 



0.0000 0000 0000 

1.3333 3333 3333 



1.0000 0000 0000 

0.3333 3333 3333 



Five-point Formula of Degree 7 

(Lobatto 5-pt Rule) 

E = -3.60 

n 

x 10- 7 f (8) (<f) 

± X. 
1 

w. 

1 



0.0000 0000 0000 

0.7111 1111 1111 



0.6546 5367 0708 

0.5444 4444 4444 



1.0000 0000 0000 

0.1000 0000 0000 



Nine-point Formula (Degree 13) 


E n = -6.16 

x 10- 16 f< 14 >(£) 

± X 

i 

w. 

1 



0.0000 0000 0000 

0.3437 6208 7210 



0.3409 8226 5911 

0.3342 3373 9816 



0.6546 5367 0708 

0.2839 7877 8048 



0.8904 0552 7513 

0.1792 6269 9553 



1.0000 0000 0000 

0.0306 4373 8977 



Seventeen-point Formula (Degree 25) 


E = -2.35 

n 

x 10- 35 f< 26) (f ) 

± X. 
1 

w i 



0.0000 0000 0000 

0.1588 8791 8790 



0.1643 3743 9165 

0.1726 1252 6200 



0.3409 8226 5911 

0.1761 6511 9745 



0.5089 0120 1269 

0.1572 3704 4683 



0.6546 5367 0708 

0.1356 0287 5607 



0.7824 1604 1120 

0.1196 2516 8318 



0.8904 0552 7513 

0.0939 8151 8229 



0.9661 9986 3094 

0.0560 0038 0783 



1.0000 0000 0000 

0.0093 3140 7070 
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Table II 


Method II 


Two-point Formula of Degree Three 

(2-point Gauss) 

E n = 7.41 

x 10- 3 f <4) (f ) 

± X. 
1 

W i 



0.5773 5026 9190 

1.0000 0000 0000 



Five-point Formula of Degree Seven 

(Kronrod) 

E = -9.00 

n 

00 

00 - 
1 

O 
T— l 

X 

± X. 
1 

w. 

i 



0.0000 0000 0000 

0.6222 2222 2222 



0.5773 5026 9190 

0.4909 0909 0909 



0.9258 2009 9773 

0.1979 7979 7980 



Eleven-point Formula of Degree 17 


E = 5.16 

n 

x 10 -22 f( 18 >(^) 

± X 

l 

w 

i 



0.0000 0000 0000 

0.3102 6774 3075 



0.3047 6216 9153 

0.2938 2188 8625 



0.5773 5026 9190 

0.2466 4959 0846 



0.7900 8322 3401 

0.1757 8249 0943 



0.9258 2009 9773 

0.0961 0769 3081 



0.9878 5368 2853 

0.0325 0446 4966 



Twenty-three-point Formula of Degree 35 

E = -1.03 

n 

x 10 _53 f (36) (<f ) 

± X. 
1 

W. 

1 



0.0000 0000 0000 

0.1551 3596 2892 



0.1544 4378 2428 

0.1530 6173 0468 



0.3047 6216 9153 

0.1469 0821 8576 



0.4469 6796 4050 

0.1368 8253 6133 



0.5773 5026 9190 

0.1233 3089 7738 



0.6926 1884 7335 

0.1067 5329 7040 



0.7900 8322 3401 

0.0878 6743 0877 



0.8679 4828 6028 

0.0677 8753 4352 



0.9258 2009 9773 

0.0481 9474 7432 



0.9650 2643 1039 

0.0306 0569 3330 



0.9878 5368 2853 

0.0155 9062 3275 



0.9978 9629 3295 

0.0054 4930 9333 







Similarly, let Q p (x) be a polynomial of degree p with zeroes at the p abscissae 
that are to be added, where 


Q P (x) =22 a j xi - 

j=o 


Then, the next formula in the sequence has n + p abscissae located at the 
zeroes of G (x) where 

n+p v 7 


and where 


G , (x) 

n+p v 7 


n + p 

•L 

j =0 



i = 0, 1 n + p 


where q = max (0, i-n) and r = min (p, i). In order to eliminate any ambiguity, 
let us assume that a Q , b 0 , and c Q are equal to one. In order to optimize the order 
or degree of the new formula, it is required that the conditions 


G (x) x k dx = 0 

n+p v ' 


-1 


k = 0, 1, . . . , p - 1 


be satisfied (Patterson, 1968). These conditions comprise a system of linear 
equations in the p unknowns { a. }p = 1 . For the two methods presented here, 
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n + p is always an odd number and for such cases the system of linear equations 


is 


n+p 

0 = y c./ (j +k + l) k = 0, 1, .... p-1. 
i^o 

j +k even 

So long as the coefficients, b. , are rational numbers, then so are the and the 
c . . Therefore, if the abscissae of the initial formula are the zeroes of a poly- 
nomial with rational coefficients, all subsequently added abscissae will also be 
the zeroes of a polynomial with rational coefficients. This fact permits the use 
of rational arithmetic in the derivation of the methods presented here in order to 
avoid the propagation of numerical errors from one stage of the process to the 
next. The only errors committed at a given stage are those arising from the 
determination of the new abscissae as the zeroes of these polynomials and from 
the computation of the corresponding weights of the quadrature formula. The 
secant method is used to obtain the zeroes of the polynomials in such a way that 
each of the zeroes is determined to lie within an interval less than 10 ~ 14 in 
length. Gaussian elimination is then used in order to determine the weights from 
the linear equations, 



= (1 - [-U j )/(j + 1) 


j = 0, 1, 


n - 1. 
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RESULTS, COMPARISONS, AND CONCLUSIONS 

The two methods presented in Tables I and II have been applied to a few 
simple quadrature problems for which the true value of the integral is known. 
The results are compared with those obtained using the following well-known 
methods: Romberg’s method, the Clenshaw- Curtis method, Patterson's method, 
and a method composed of Gauss formulae. The composite Gauss method con- 
sists of Gauss formulae of degrees one, three, seven, and fifteen. The total 
number of function evaluations required at each stage are therefore one, three, 
nine, and twenty- three respectively. 

The following set of problems are considered: 


/•i 

8 dx 

*Li x 2 + 2x+5 


= rr, 


( 1 ) 


'l 

x 1/2 dx = 2/3, 
J o 


( 2 ) 


'i 
« . 


x 3/2 dx = 2/5, 


(3) 



l 


(4) 


and 


I dx 
h I + x 2 


tt/ 2 . 


(5) 
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The results of applying the above methods to this set of problems are 
summarized in Table III. 

The first three methods shown in the table use the function values at the 
endpoints of the interval of integration and the last three methods do not. For 
this reason, let us examine these first three methods as a group. In problems 
1, 4, and 5, the method which results in the smallest error is Method I, followed 
by the Clenshaw- Curtis method, and Romberg’s method. In examples 2 and 3, 
which possess singular derivatives at the origin (first and second derivative, 
respectively), the Clenshaw- Curtis method produces the best approximation, 
followed by Method I and the Romberg method in that order. 

For the last three methods in the table, the comparison is more difficult 
because the different methods produce results for different values of n, the total 
number of function evaluations. In any case, if we consider the number of sig- 
nificant figures gained per function evaluation as a measure of the efficiency 
of a method, some conclusions can be drawn. In each of the five problems, the 
Patterson method and Method II provide results that are almost identical using 
the above measure of efficiency. The results using the composite Gauss method 
are somewhat inferior. 

Comparison of all six methods tends to indicate that the last three methods 
are more efficient, especially on problems 2 and 3. In all cases, the composite 
Gauss method performed as well or slightly better than the Clenshaw- Curtis 
method. For each problem, the best methods are Patterson’s method and 
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Table m 


Errors Obtained Using Certain Quadrature Methods 


Problem 

n 

Romberg's 

n 

Clenshaw- Curtis 

n 

Method I from 

n 

Patterson's 

n 

Composite Gauss 

n 

Method II from 

No. 

Method 

Method 

Table I 

Method 

Method (1-3-7-15) 

Table II 


3 

8.26 x 10~ 3 

3 

8.26 x 10‘ 3 

3 

8.26 x 10' 3 

i 

5.84 x lO" 2 

i 

5.84 x lO* 2 

2 

5.94 x 10‘ 3 

1 

5 

5.25 x 10 ' 4 

5 

2.36 x 10~ 4 

5 

2.40 x 10‘ 5 

3 

5.25 x 10' 4 

3 

5.25 x 10~ 4 

5 

6.04 x 10‘ 6 


9 

6.87 x 10 " 6 

9 

3.31 x 10' 9 

9 

6.91 x 10' 10 


1.84 x 10‘ 8 

9 

2.66 x 10' 9 

11 

6.32 x 10" 13 


17 

1.17 x 10“ 8 

17 

- 

17 

- 

l 

- 

23 

- 

23 

- 


3 

2.86 x 10 ' 2 

3 

2.86 x 10‘ 2 

3 

2.86 x 10- 2 

B 

4.04 x lO* 2 

1 

4.04 x 10 " 2 

2 

7.22 x 10' 3 


5 

8.91 x 10' 3 

5 

2.07 x 10- 3 

5 

4.53 x 10 - 3 

m 

2.51 x 10' 3 

3 

2.51 x 10‘ 3 

5 

3.41 x 10 - 4 

2 

9 

3.06 x 10 " 3 

9 

2.24 x 10- 4 

9 

7.61 x 10 * 4 


1.42 x lO' 4 

9 

2.46 x 10 ~ 4 

11 

2.48 x 10 ~ 5 


17 

1.07 x 10 " 3 

17 

2.69 x 10- s 

17 

1.27 x 10' 4 

B 

7.05 x 10 ‘ 6 

23 

2.77 x lO’ 5 

23 

2.13 x 10‘ 6 


3 

2.37 x 10~ 3 

3 

2.37 X 10- 3 

3 

2.37 x 10~ 3 

i 

4.64 x 10 - 2 

1 

4.64 X 10 * 2 

2 

1.22 x 10' 3 


5 

3.03 x 10~ 4 

5 

1.25 x 10 * s 

5 

9.21 x lO* 5 

3 

1.88 x 10~ 4 

3 

1.88 x 10- 4 

5 

2.76 x lO" 6 

o 

9 

4.96 x 10 ~ S 

9 

7.79 x 10 - 7 

9 

4.53 x lO’ 6 

7 

1.03 x 10 “ 6 

9 

3.60 x 10* 6 

11 

7.71 x 10'® 


17 

8.62 x 10' 6 

17 

2.73 x 10 ' 8 

17 

2.27 x 10' 7 

15 

1.78 x 10 -9 

23 

9.29 x lO' 8 

23 

1.63 x lO -9 


3 

1.30 x 10- 3 

3 

1.30 x 10' 3 

3 

1.30 x lO" 3 

1 

2.65 x 10~ 2 

1 

2.65 x 10 - 2 

2 

8.39 x 10' 4 

4 

5 

2.74 x 10* 5 

5 

9.93 x 10' 6 

5 

9.68 x 10' 7 

3 

2.55 x lO" 5 

3 

2.55 x 10 " 5 

5 

2.18 x 10' 7 


9 

2.97 x 10‘ 7 

9 

6.40 x 10' 10 

9 

6.44 x lO’ 12 

7 

2.39 x 10 ~ 10 

9 

1.99 x 10* u 

11 

- 


17 

1.36 x 10~ 9 

17 

- 

17 

- 

15 

- 

23 

- 

23 

- 


3 

9.59 x 10" 2 

3 

9.59 x 10~ 2 

3 

9.59 x 10- 2 

1 

4.29 x lO' 1 

1 

4.29 x 10 - 1 

2 

7.08 x 10- 2 


5 

1.08 x 10- 2 

5 

6.98 x 10' 3 

5 

2.54 x 10- 3 

3 

1.25 x 10 - 2 

3 

1.25 x lO' 2 

5 

9.99 x 10 ' 4 


9 

4.38 x 10- 4 

9 

9.77 x 10' 6 

9 

1.06 x 10' 6 

7 

3.35 x 10 - 5 

9 

1.11 x 10 -s 

11 

1.37 x 10’ 7 


17 

5.17 x 10‘ 6 

17 

5.30 x 10 ' 10 

17 

4.88 x 10-n 

15 

1.87 x lO' 10 

23 

8.46 x 10- 12 

23 

- 


— Indicates an error less than 1 x 10-' 3 n = Total number of function evaluations. 



































Method II, each giving about the same over-all efficiency. Following these for 
examples 2 and 3 are the composite Gauss method, the Clenshaw- Curtis method, 
Method I, and Romberg's method, in that order. For examples 1, 4, and 5, 
Method I is the third most efficient method, followed by the composite Gauss, 
the Clenshaw- Curtis method, and Romberg's method, in that order. 

These preliminary results seem to indicate that the methods of the Patter- 
son type are certainly worthy of further development and study. 
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